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S, Rangaswamy 


ABSTRACT 


Two techniques for recovering atmospheric refractivity profiles from 
simulated satellite-to- satellite trackiilg data are documented and the 
relative degradation in accuracy due to the calculation methods shown 
(modelling errors not included). Examples are given using the geo- 
metric configuration of the ATS-6/ NIMBUS-6 Tracking Experiment. 
Both satellites are in circular orbit (NIMBUS in polar orbit at a height 
of 1000 kilometers and ATS in a geostationary equatorial orbit). The 
vmderlying refractivity model for the lower atmosphere has the spheri- 
cally symmetric form N = exp P(s) where P(s) is a polynomial in the 
normalized height s. For the particular simulation used in the analysis, 
the Herglotz-Wiechert technique recovered values which were 0.4% and 
40% different from the input values (calculated from the parameters of 
the assumed model atmosphere) at the surface and at a height of 33 ki- 
lometers, respectively. Using the same input data, the model fitting 
technique recovered refractivity values 0.05% and 1% different from 
the input values at the surface and at a height of 50 kilometers, 
respectively. 

It is also shown that if ionospheric and water vapor effects can be 
properly modelled or effectively removed from the data (e.g. , by sepa- 
rate measurement such as a radiometer or by regression techniques), 
pressure and temperature distributions can be obtained from the dry 
refractivity by numerical integration, assuming hydrostatic equilibrium 
and the perfect gas law. This is demonstrated for radiosonde data 
taken at Dulles airport, where pressure is recovered to less than 0.3% 
up to 20 kilometers. The temperature recovery for the particular ex- 
ample was 5% at the surface but improved with height. 


iii 


CONTENTS 


Page 

1.0 INTRODUCTION 1 

2.0 TRACKING MEASUREMENTS 3 

3.0 LOWER ATMOSPHERIC REFRACTIVITY PROFILE 

AND MODEL * 4 

4.0 EFFECT OF LOWER ATMOSPHERE ON RANGE, 

RANGE RATE, AND BENDING ANGLE . . ...... 5 

5.0 IONOSPHERIC MODEL AND EFFECT OF IONOSPHERE 

ON TRACKING DATA .13 

6.0 ITERATIVE MODEL FITTING TECHNIQUE. ..... . 14 

7.0 THE HERGLOTZ-WEICHERT METHOD ................... 22 

8.0 RECOVERY OF PRESSURE AND TEMPERATURE ............ 28 

9.0 SUMMARY 33 

10. CONCLUSIONS 34 

REFERENCES. . . . . . ... 35 

APPENDIX A ........ . . .... . . . . . . A-1 

APPENDIX B . . .... ... . . B-1 

APPENDIX C ... . . C-1 


PRECEDING PAGE BLANK NOT PIU® 


V 



ILLUSTRATIONS 


Figure Page 

1 ATS-6/NIMBUS-6 Geometry 2 

2 Tracking Measurements 3 

3 Lower Atmospheric Refractivity Profile 6 

4 Residuals from a 7th Degree Least Squares Polynomial Fit 

to the Logarithm of Refractivity 7 

5 RMS of Least Squares Polynomial Fits to the 

Logarithm of Refractivity. 8 

6 Lower Atmospheric Effects for Typical Refractivity Profile . . . . 9 

7 Difference Between Atmospheric Range Error Along a 

Straight Line and Along the Ray Path of the Signal 12 

8 Comparison of Ionospheric Electron Density Models 16 

9 Ionospheric Effect for Typical Ionospheric Density 17 

10 Occultation Geometry for Iterative Model- Fitting Technique . ... 20 

11 Flow Chart of Iterative Model- Fitting Algorithm 21 

12 Occultation Geometry for Herglotz-Wiechert Technique ... ..... 24 





VI 



TABLES 


Table Page 

1 Comparison of Various Ionospheric Models Fit to Composite 

Electron Density Profile in Figure 8 15 

2 Recovered Lower Atmospheric Refractivity by 

Iterative Model Fitting 23 

3 Results of hiversion by H-W Algorithm for the 

Three Parameter Profile 27 

4 Recovery of Pressure and Temperature , . 32 


vli 


RECOVERY OF ATMOSPHERIC REFRACTIVITY PROFILES 
FROM SIMULATED SATE LLITE-TO-SATE LLITE TRACKING DATA 


1.0 INTRODUCTION 

The purpose of this report is two-fold: (1) to document two techniques for re- 
covering lower atmospheric refractivity profiles^ from simulated satellite-to- 
satellite tracking data, and; (2) to show how pressvire and temperature can be 
calculated from the recovered refractivity profile provided the effects of the 
ionosphere and water vapor are properly modelled or effectively removed from 
the data. Also shown are the effects of the atmosphere upon the one-way range 
and range rate between the two satellites. 

Both the geometry of satellite-to-satellite tracking and precise range and range 
rate measurements expected in the ATS-6/NIMBUS-6 Tracking and Orbit Deter- 
mination Experiment [ 1] have provided the motivation for the present analysis. 

The geometry used for the simulation can be seen in Figure 1 where NIMBUS is 
in a circular polar orbit of 1000 kilometers and ATS is in a geostationary orbit 
in the equatorial plane. ^ 

The tracking signal travels a 4-way path: from the grotind station to ATS, from 
ATS to NIMBUS, from NIMBUS to ATS, and ATS back to grovuid. 

Radio occidtation occurs whenever the signal between the two satellites passes 
through the atmosphere. As NIMBUS emerges from or enters into radio shadow 
behind the earth, the phase of the tracking signal between the two satellites is 
affected by the presence of the earth's atmosphere. Given the ephemeris of 
each satellite and the tracking data (range and/or range rate), it is possible to 
obtain a vertical refractivity profile of the earth's lower atmosphere above the 
point where the ray path between the satellites comes closest to the earth. 


^For this analysis the lower atmosphere is defined to be the neutral atmosphere from the surface of the 
earth to an altitude of 100 kilometers. The ionosphere is defined as that part of the atmosphere containing 
charged particles in sufficient density to affect the propagation of S-band signals (2 GHz) and extending 
from an altitude of 100 kilometers to 1000 kilometers. 

^NlMBUS-6 was successfully launched 12 June 1975 and is in a near-circular orbit (eccentricity of 'V/O.OOOS) 
at an altitude of 1000 kilometers and inclination of approximately 100°. ATS-6 is also in a near-circular 
orbit (eccentricity of 'vO.OOOS) with an inclination of 0.9° and semi-major axis of about 41870 kilometers. 
The orbits of ATS and NIMBUS are assumed to be circular for simplicity in analysis. 
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z (NORTH) 



NOTE: FIGURE DOES NOT DEPICT NOMINAL OR ACTUAL ORBITS. 
CIRCULAR ORBITS ASSUMED 

Figure 1, ATS-6/NIMBUS-6 Geometry 


The mathematical techniques of profile inversion are well- documented [2] and 
planetary atmospheres have been studied using microwave radio occultation 
methods for Mars ([3] , [4] , [5]) and Venus ([6] , [7] , [ 8] ). In [9] it is pointed 
out that the ATS-6/NIMBUS-6 Tracking Experiment provides a ready-made op- 
portunity for measuring atmospheric profiles of the earth. Marini [9] assumes 
an exponential model for the lower atmosphere and a parabolic model for the 
ionospheric density and gives the errors in range and range rate between the two 
satellites for a 2 GHz carrier signal. The calculations are made for straight 
line propagation and show measurable atmospheric effects of 150 meters in range 
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and 50 meters/ second in range rate and ionospheric effects of 15 meters in 
range and 30 centimeters/ second in range rate. 


Both the Herglotz-Wiechert method and the iterative model-fitting approach 
which are described in this analysis can be found in the literature [ 2] . They 
have been adapted for use here in order to meet the geometrical constraints 
imposed by the ATS-6/NIMBUS-6 Tracking and Orbit Determination Experi- 
ment [ 1] . 


2.0 TRACKING MEASUREMENTS 

For the ATS-6/NIMBUS-6 Tracking Experiment [1] , the signal travels a 4-way 
path as shown in Figure 2 which is taken from [ 1] . The range measurement 
consists of the time in seconds between a zero crossing of the ground trans- 
mitted signal (6 GHz carrier phase modulated by a 100 KHz sidetone) and a zero 


ATS-F 



ROSMAN, N.C. 
tracking STATION 

Figure 2. Tracking Measurements (taken from [1] ) 



crossing of the ground received signal (4 GHz carrier having two subcarriers 
each phase modulated by the 100 KHz sidetone) after the signal has travelled 
from groimd to ATS, to NIMBUS, back to ATS, and then to ground ({ 10] ). The 
range rate measimement consists of the time in seconds required to count a 
fixed number of cycles of the 4-way Doppler frequency plus a fixed offset fre- 
quency. (De struct Mode, [10] ,) 

For purposes of this report the "observables” are the one-way range between 
NIMBUS (at the time the signal leaves there) and ATS (at the time the signal 
arrives there) alorit the propagation path of the signal and the rate of change of 
this quantity with time or range rate along the propagation path. 


3.0 LOWER ATMOSPHERIC REFRACTIVITY PROFILE AND MODEL 
Refractivity at radio frequencies is given by ([11] , pg. 7) 

N = 77.6 (P/T) + (3.73) 10® (e/T2) (3.1) 

where P is pressure in millibars, T is temperature in degrees Kelvin, and e is 
the partial pressure of the water vapor in millibars given by ([12] , pg. 343) 

7.5 (T - 273. IS) 

e = (R^ / 100) (6. 11) 1025^' 3 + (T- 273.15) ^3 2) 

where R^ is the relative humidity in percent. ® 

Refractivity N is related to the index of refraction n^ by the following 

N = (n - 1) 10® (3.3) 


The first term on the right-hand side of (3.1) is the dry refractivity N^ and the 
second term on the right-hand side the wet refractivity N^. 

the dewpoint tempera^jre (°K) is given instead oi‘ the relative humidity, e can be calculated from 
(3.2) by setting Rj, equal to 100 and equal to T. 

^The group and phase index of refraction are the same since the lower atmosphere is assumed to be non- 
dispersive. 
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An example of a ref r activity profile calculated from tables of pressime, tempera- 
ture, and relative humidity versus height for a standard atmosphere of 30° North 
July [13] , using (3.1) can be seen in Figure 3. Both the dry and the dry-plus-wet 
profiles are shown. 

It can be seen that log^N is almost linear in the height h. Bean and Thayer [14] 
have found the average variation of N with height to be nearly exponential. 

Hence, it seems reasonable to model the lower atmospheric refractivity with 


N = exp P(s) 


(3.4) 


where P(s) is a polynomial in the normalized height s = (.01) h - 1 (h is height 
in kilometers). 

Figure 4 shows the residuals from a seventh degree least squares polynomial 
fit to loggN where N is the dry-plus-wet refractivity shown in Figure 3. 

Figure 5 shows the increased goodness-of-flt with higher degree polynomial fits 
to loggN. Both the dry and the dry-plus-wet refractivities are compared and it 
is seen that there is little difference between them as the degree is increased. 


4.0 EFFECT OF LOWER ATMOSPHERE ON RANGE, RANGE RATE, 

AND BENDING ANGLE 

A radio wave passing through the earth's atmosphere is bent and, in addition, 
delayed along its path. The first effect is due to the varying index of refraction 
across the signal path, while the second effect is due to the fact that the index of 
refraction for the lower atmosphere is greater than unity. Both these effects 
combine to cause the apparent range between a transmitter and receiver to be 
greater than the line-of-sight range. 

Figure 6 shows the lower atmospheric effects upon range, range rate, and bending 
angle for a realistic refractivity profile of the form indicated in (3.4) where P(s) 
is a second degree polynomial. The coefficients of P(s) were determined using 
linear least squares, fitting the pol 3 momial to the logarithm of the calculated 
values of refractivity. The refractivity was calculated from radiosonde data 
taken at Dulles airport, Virginia, in 1967, using (3.1) and the procedure outlined 
in [15] . The geometric configuration of ATS and NIMBUS is indicated in the 
figure where ATS is near the orbital plane of NIMBUS. 
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Figure 3. Lower Atmospheric itefractivity Profile 
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Figure 4. Residuals from a 7th Degree Least Squares Polynomial Fit to the Logarithm of Refractivity 



F,ms X W WHERE IS IN UNITS OF LOG N (N = REFRACTIVITY) 


LEGEND: - 


DRY REFRACTIVITY 

DRY PLUS WET REFRACTIVITY 



NOTE: THE REFRACTIVITY PROFILES WHICH WERE FIT ARE 
SHOWN IN FIGURES. 

REFRACTIVITY N WAS CALCULATED FROM TABLES OF 
TEMPERATURE, PRESSURE, AND RELATIVE HUMIDITY FOR 
A STANDARD ATMOSPHERE OF 30° NORTH JULY ([13]). 
THE FORMULA USED FOR DRY-PLUS-WET REFRACTIVITY IS 

N = No + Nyv 

WHERE THE DRY REFRACTIVITY Np IS GIVEN BY: 

Nd = 77.6 (P/T) AND WET BY Nw = 3.73 X (e/T^), 

WHERE P AND e ARE IN MILLIBARS AND T IN °KELVIN. 




DEGREE OF LEAST SQUARES PDLYNOMIAL 


Figure 5. KMS of Least Squares Polynomial Fits to the Logarithm of Refractivity 
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Figure 6. Lower Atmospheric Effects for T 3 ?pical Refractivity Profile 


The range error due to the lower atmosphere (see Figure 6) is 



where 


(4.1) 


- Ri + Ro + 2 




ds 


or 




R, = Ri + R 2 + 2 


Rn - Ra 


ndr 


|/l - (noro/nr)^ 


(4.2) 


and Ri and Rj are the linear segments taken from the "top^' of the lower atmo- 
sphere r^ (at an altitude of 100 kilometers) to the radius of ATS and the radius 
of NIMBUS, n is the index of refraction at the geocentric radius r, R^ and R^ 
are the position vectors of ATS and NIMBUS, R^ is the line-of-sight range. 

The bending angle T [16] is given by 


- = -r 


cot d (dn/n) 


or 


r- m 


r "" 2riQ T0 10 ^ 


(dN/dr) dr 


n^ r y'l - (ii(,r(,/nr)‘ 


(4.3) 
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using (3.3) and Snell's law for a spherically stratified medium 


nr cos 6 = iig Tq cos 00 (4.4) 

where d is the local elevation angle of the ray at radius r and n is the index of 
refraction at the radius r. 

The error in range rate due to the lower atmosphere is 

AR = - R^ (4.5) 


where is the range rate along the propagation path which is obtained by sum- 
ming the dot products of the velocity vectors of ATS and NIMBUS with unit vec- 
tors along the propagation path at each satellite (using the convention that an 
increasing range with time gives rise to a positive range rate). The unit vectors 
are obtained through ray trace procedures and there is a one-to-one correspon- 
dence between these vectors and the rai^e • R^ is the range rate along the 
line-of« sight range R^. 

The quantities AR, AR, and r are obtained as functions of time by linking them 
to the ephemeris of both satellites. 

From the figure the atmospheric range error is more than 2 kilometers and the 
range rate error more than 100 meters/ second at the time of occultation. This 
is largely due to the bending angle of approximately 2° at the surface. 

Figure 7 shows the difference between the atmospheric range error taken along 
a curved path and the range error taken along a straight line path. Also shown 
is the difference between these two quantities which is 38 meters at the time of 
occultation. 

The assumption of a straight ray (neglecting bending) simplifies the analysis and 
is valid for tenuous atmospheres like that of Mars. However, Figures 6 and 7 
show that in the case of the earth's atmosphere, refractive bending of the ray 
cannot be neglected. 
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Figure 7. Difference Between Atmospheric Range Error Along a Straight Line 

and Along the Ray Patli of the Signal 
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5.0 IONOSPHERIC MODEL AND EFFECT OF IONOSPHERE 
ON TRACKING DATA 

Although the main purpose of this analysis is the recovery of lower atmospheric 
refractivity profiles, the ionosphere caimot be neglected at S-band frequencies 
(2 GHz). Its effect, however, is not as great as the lower atmosphere upon one- 
way satellite-to-satellite range and range rate [9] . The manner in which the 
ionosphere enters into the occultation calculations is in the computation of the 
path length of the signal through that region (100 to 1000 kilometers). See equa- 
tion (B.9) of Appendix B. But for rays passing close to the earth, the change in 
the contribution of the ionosphere to the total path length should be rather small. 
For this reason and the fact that the ionospheric errors at S-Band are much less 
than the lower atmospheric errors (in the tracking data), a simple model for the 
ionospheric density is desirable. The ionosphere is dispersive and therefore the 
distinction must be made between group and phase velocities and group and 
phase indices of refraction. The group index of refraction is related to the group 
refractivity (equation 3.3) by 


Ng = (tig - 1) 10^ (5.1) 

and Ng is given by [17] 


40.3N^(h) 

Ng(h) = —7— 10" (5-2) 


where 

N^ - electron density (electrons/meter^) 
f = frequency (Hz) 
h = height 

With the objective of simplicity in mind a Chapman and a modified Chapman 
function [18] were fit to a composite profile from Alouette-I and Wallops Island 
ionograms [19].^ The results can be seen in Figure 8. The modified Ch^man 
function yields a much better fit than the Chapman function up to an altitude of 
about 700 kilometers. ^ ^ 
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Table I shows a comparison of various ionospheric models fit to the composite 
profile in Figure 8. It can be seen that the exponential polynomial gives the best 
fit (Columns 3 and 4). 

Figure 9 shows the effect of the ionosphere for a typical electron density 
profile. 

The effect of the ionosphere in the recovery of lower atmospheric refractivity 
profiles is included in the next section. 


6.0 ITERATIVE MODEL FITTING TECHNIQUE 

The problem here is to determine the parameters of an assumed refractivity 
model for the lower atmosphere of the form (3.4) given an ephemeris of the two 
satellites and the one-way range and/or range rate between them for a number 
of data points. The technique uses the generalized method of nonlinear weighted 
least squares with apriori parameter vector x and apriori covariance matrix 
(Appendix A) to find the (m x 1) parameter vector 

= (ag ai 32 .. . (6.1) 

of the poljniomial P 

P = s'^'x = s + . . . + (6*2) 

^ The Chapman profile is given by 

N(h) = Nmaxexp ^ (1 -z-e-^) 

where 

N(h) = electron density at an altitude h(e/m^) 

Nmax ~ maximum electron density (e/m^) 

Z = (h-hniax)/H 

hmax ^ height at which electron density is maximum 
H = scale height 

The modified Chapman profile is given by 

N(h) = Nmaxexp \ 

A = 1 - exp 
Z' = (h-h„,ax)/Ho 

the parameters of the model being Nmax'bmax'Q' Hq . 
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Table I 


Comparison of Various Ionospheric Models Fit to Composite 
Electron Density Profile in Figure 8 


DAY 323, 1970 LMT 1549 SOL. ZEN. ANGLE 80.5° 


Altitude 

(km) 

Electron Density (10 Electrons/m^) 

Composite 

Profile 

Polynomial 
Degree 7 

Polynomial 
Degree 10 

Ch^man 

Modified 

Ch^man 

200 

0.41 

0.40 

0.40 

0.50 

0.39 

250 

1.29 

1.29 

1.30 

0.93 

1.29 

275 

1.43 

1.46 

1.43 

1.02 

1.29 

300 

1.34 

1.32 

1.33 

1.02 

1.15 

325 

1.13 

1.11 

1.12 

0.95 

0.99 

350 

0.91 

0.90 

0.91 

0.86 

0.84 

400 

0.58 

0.58 

0.58 

0.65 

0.59 

450 

0.38 

0.39 

0.39 

0.46 

0.42 

500 

0.27 

0.27 

0.27 

0.32 

0.30 

550 

0.21 

0.21 

0.21 

0.22 

0.22 

600 

0.16 

0.16 

0.16 

0.15 

0.16 

650 

0.13 

0.13 

0.13 

0.10 

0.12 

700 

0.10 

0.10 

0.10 

0.07 

0.09 

750 

0.08 

0.08 

0.08 



800 

0.07 

0.07 

0.07 



850 

0.06 

0.06 

0.06 



900 

0.05 

0.05 

0.05 



950 

0.05 

0.05 

0.05 



1000 

0.04 

0.04 

0.04 
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Figure 9. Ionospheric Effect for Typical Ionospheric Density 


in the model (3.4) 


N - exp P(s) 


(6.3) 


which minimizes the quadratic form 

Q = [z - f (x)]^ [z - f (X)] + (X - x)^ Djci'(x - x) (6.4) 


The vector s 


sT = (1 s s2 ... s'"-!) (6.5) 

is the vector of normalized heights s = (0.01) s - 1 (h is the height in kilometers). 
The vector x 


X - (aoaia.^ ... a„_^) 


( 6 . 6 ) 


is the apriori vector of assumed starting values for the coefi'icients aj in (6.2) 
with associated covariance matrix Dji 

The vector z 

z-r - (zjzj ... z„) (6.7) 


is the observation vector (range and/or range rate values) with covariance 
matrix (taken to be the identity matrix in this report, i.e., least squares with 
equal weights). 


The vector f 



is the vector of calculated quantities (range and/or range rate) which is anon- 
linear function, of the vector x and is obtained by ray tracing the radio signal 
connecting the positions of ATS and NIMBUS through the assumed model (6.3) 
with the initial estimate of the parameters x in (6,6). By comparing the com- 
puted values of the range or range rate f in (6.8) with the observed values z in 
(6.7), it is possible to correct the parameter set iteratively. 

As shown in Appendix A, using a modified Newton-Raphson iteration scheme, 
the solution vector x* is 


X* = x(>">+ [bTD-^B + |bTd-^ [iz- f (xW)] + (x - x(">))| (6.9) 

where B is an (n x m) matrix of partial derivatives (rank m < n) 


= B(x) (6.10) 

and X is the estimate of the parameter vector x at the m^ iteration. 

The variance of the refractivity at the normalized height s is given by 



cr^ - exp (2 X*) [B^D“'B + D~ s 


( 6 . 11 ) 


The standard deviation of N is the square root of (6,11). 

It shovild be noted that the location of the specific ray connecting the two satel- 
lites is in itself an iterative procedure. Starting with a given radius of closest 
approach rp (Figure 10), a ray is traced through the given atmosphere using 
Snell's law for a spherically stratified medium (4,4). At the assumed limit of 
the atmosphere (at a height of 100 kilometers for this analysis), the ray tracing 
is stopped and the ray asymptote is extended until it intersects the circle with 
radius equal to the radius vector of the satellite. The angle subtended by the 
two points of intersection at the center of the eax^ is compared with the angle 
subtended by the locations of ATS and NIMB US at the center of the earth. If the 
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Figure 10. Occultation Geometry for Iterative Model-Fitting Technique 


two angles do not agree to the required precision (10"® radians), the radius of 
closest approach is changed and another ray is traced. A modified Newton- 
Raphson iterative procedure is used. 

After die ray connecting the two satellites is found, the range along that ray, as 
well as the partial derivatives of the range with respect to the parameters of the 
model atmosphere are determined, (6,10) (in the case where the observations are 
ranges). The partial derivatives are found numerically. This procedure is re- 
peated for all the given pairs of positions of the two satellites. Computed values 
of range f are compared with observed values z. If they do not agree to the re- 
quired precision, the residuals along with the partial derivatives are used to 
estimate corrections to the parameter set x by using the least squares criterion 
(Appendix A). This process is repeated iteratively imtil the change in the parame- 
ter set is within a specified magnitude. The solution is then given by (6.9). Fig- 
ure 11 gives the flowchart of the computation. 

Results of Iterative Model Fitting 

Simulated range data along with the sequential position and velocity coordinates 
of ATS and NIMBUS were used as input for the iterative model fitting technique 
described above. In all cases investigated the range data included both lower 
atmospheric and ionospheric effects. The input lower atmospheric refractivity 
profile was a least squares fit of N - exp [ P(s)J to values of refractivity versus 
height calculated from radiosonde data taken at Dulles airport, Virginia, during 
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Figure 11. Flow Chart of Iterative Model- Fitting Algorithm 

(Range Data) 
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1967 where P(s) is a second degree polynomial. The input ionospheric electron 
density profile is a least squares fit of exp [Q(s)] to the composite profile indi- 
cated in Figure 8 where Q(s) is a seventh degree polynomial (see Table I). The 
results of the model fitting technique can be seen in Table II. Colvunn 1 is the 
height in kilometers above the surface. Column 2 is the input lower atmospheric 
refractivity. In all cases the lower atmospheric refractivity profile which is 
’’solved for” is a function of the form exp (P(s)] where P(s) is a 2nd degree 
polynomial. 

The third column is the recovered refractivity where the input ionospheric elec- 
tron density function (8 parameter model of form exp [Q(s)]) has been ’’modelled” 
into the calcxilations. The error in recovery is on the order of 0,05% at the sur- 
face and less than 1% at 50 kilometers. 

The fourth column shows the results of neglecting the ionosphere entirely 
(neither ’’modelled” nor ’’solved for”). The errors in recovery are large; about 
7% at the surface and about 375% at around 45 kilometers height. This clearly 
indicates that the ionosphere cannot be neglected at S-band frequencies. 

The fifth column shows the resvilts of modelling the Chapman function in the 
calculations (the parameters for this model were obtained by fitting in the least 
sense the Chapman function to the composite profile in Figure 8). The error in 
recovery is about 1.3% at the surface and 25% at 40 kilometers. 

The sixth column is the recovered refractivity where the ionospheric electron 
density has been ’’solved for” and has the Chapman form. The maximum error 
is about 12.5% at 40 kilometers. 

Finally, in the last column are shown the recovered refractivity values when the 
parameters of the modified Ch^man function are ’’solved for.” Recovery in 
this case is excellent and is seen to be approximately 0.3% at 15 kilometers 
(maximum). 

The above demonstrates: (1) the ionosphere cannot be neglected as S-band fre- 
quencies (colvunn 4), and; (2) a reasonably simple model for the ionospheric 
electron density such as the modified Chapman should be sufficiently accurate to 
recover lower atmospheric refractivity profiles (column 7). 

7.0 THE HERGLOTZ-WIECHERT METHOD® 

This technique was originally developed for the pui^jose of inverting seismic 
data in order to obtain velocity depth profiles of the earth. The analysis makes 
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Table n 


Recovered Lower Atmospheric Refractivity By Iterative Model Fitting 


Height 

(Kilometers) 

Input Lower 

Atmospheric 

Refractivity 

Recovered Lower Atmospheric Refractivity 

Lower Atmospheric 
Refractivity Only 
"Solved For." Eight 
Parameter Ionospheric 
Profile "Modelled." 

Lower Atmospheric 
Refractivity "Solved 
For." No Ionosphere 
"Modelled" or 
"Solved For." 

Lower Atmospheric 
Refractivity "Solved 
For." Chapman 

Function "Modelled." 



Lower Atmospheric 
Refractivity "Solved 
For." Chapman 
Function "Solved For." 

Lower Atmospheric 
Refractivity "Solved 
For." Modified 
Chapman Function 
"Solved For." 

0 

375.2 

375.4 

348.7 

370.4 

372.9 

375.3 

5 

175.3 

175.3 

195.8 

176.2 

175.9 

175.3 

10 

81.9 

81.9 

IIQ.O 

83.8 

83.0 

81.9 

15 

38.3 

38.3 

61.8 

39.9 

39.1 

38.2 

20 

17.9 

17.9 

34.7 

18.9 

18.4 

17.8 

25 

8.3 

8.3 

19.5 

9.0 

8.7 

8.3 

30 

3.9 

3.9 

10.9 

4.3 

4.1 

3.9 

35 

1.8 

1.8 

6.1 

2.0 

1.9 

1.8 

40 

0.8 

0.8 

3.4 

1.0 

0.9 

0.8 

45 

0.4 

0.4 

1.9 

0.5 

0.4 

0.4 

50 

0.2 

0.2 

0.2 

0.2 

0.2 

0.2 


Note 1: In all cases the satellite-to- satellite range data includes both lower at- 
mospheric and ionospheric effects. The input lower atmospheric refrac- 
tivity profile is a least squares fit of exp [P(s)] to values of refractivity 
versus height calculated from radiosonde data taken at Dulles Airport, 
Virginia in 1967 where P(s) is a second degree polynomial. The input 
ionospheric electron density profile is a least squares fit of exp [Q(s)] 
to the composite profile in Figure 8 (Q(s) is a 7th degree polynomial, 
see Table I). 

Note 2: la all cases above the lower atmospheric refractivity profile which is 
"solved for" (in the least squares sense) is a 3 parameter model of the 
form N = exp [ P (s)] . 

Note 3: In the above the term "modelled" (e.g. Chapman function "modelled") 

means the least squares fit of that function to the composite profile (in- 
dicated in Figure 8 and Table I) was fixed in the calculations. 
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use of the Abel transform for the integral describing the phase change along a 
bent ray path. This method has been adapted for radio occultation data by 
Phinney and Anderson [ 20] whose derivation is given in Appendix C. The re- 
fractivity index n, corresponding to a radius r is given by 



r /3 

L\ 

n ( r ) exp 

(1/tt) I oosh"^ (p/t 7) f- 

m 



where = P(/3) - 2r^cos i„, i^ is the angle of incidence at the outer limit of 
the atmosphere, and P(/o) is the phase path length (equation C.2). See Figure 12. 



Figure 12. Occultation Geometry for Herglotz-Wiechert Technique 


®This section is self-contained. The notation used is different from that used elsewhere in the report. 
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Computation of Impact Parameter 


It can be sho^ [20] that the angle between the satellite (NIMBUS) trajectory 
and the ray path ’Z' (Figure 12) is given by 


cos yjj - (D^ + Dg)/V 


(7.2) 


where is the contribution to the rate of phase path change due to the motion 
of the satellite and n, the contribution due to the atmospheric refraction effect, 
V is the projection of the velocity of NIMBUS on the plane containing NIMBUS, 
ATS, and the center of the earth. Similarly, the angle between the trajectory of 
NIMBUS and the line of sight between the two satellites is given by 

cos = Dg/V (7.3) 

The range rate between the two satellites can be shown to be the sum of the pro- 
jection of the velocity vector of each satellite along the ray path. 

When there is no atmosphere the true range rate is given by 

^AL (7.4) 


where V^l is the projection of the velocity vector of NIMB US along the line of 
sight and V., is the projection of the velocity vector of ATS along the line of 
sight. From (7.3) 


= V/V 


Combining (7.4) and (7.5) we get 


0, = cos-1 [(Rt“ Val)/v] (7.6) 
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Similarly we can show 



where R is the measured range rate through the atmosphere and Var is the 
projection of the velocity vector of ATS along the ray path. 

From Figure 12 we see that the impact parameter P is given by 


p - ]Rj^1 sin/S (7.8) 

where I Rfj l is the radial distance of NIMBUS and 


^ ^ I/; + <0N A (7.9) 

Hence, by knowing 0, , and the trajectories of NIMBUS and ATS it is possible 

to compute the impact parameter. 

However, the computation of V' requires a loiow ledge of the ray path (equation 7.7) 
But to start with, the ray path which is uniquely defined by the impact parameter, 
is not known. Hence it is necessary to undertake an iterative procedure to deter- 
mine the impact parameter. 

As a first approximation it is assumed that Var - Val . and 0, and p are 
computed. Now 


slnl„ - p/t^ /7.10) 

By knowing i^ it is possible to compute the direction cosines of the ray at 
ATS and hence Var The new value of Va^ is then input to determine a new value 
of p. This process is continued until convergence is reached. 

Computation of Refractivity 

The evaluation of refractivity from equation (7.1) is done in the following steps: 
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(1) The input data consisting of position and velocity of the two satellites 
as a function of time, the range and range rate along the propagation 
path, are read in. 

(2) For each data point the impact parameter is computed as described in 
the above. The range residual (equation 7.1) is also computed for 
each data point by subtracting the linear segments of the ray from the 
top of the atmosphere to each satellite from the range, hi addition the 
term 2 r^ cos i^^, is also subtracted from the range. 

(3) The logarithm of Pg is then fitted as a polynomial in p , normalized be- 
tween the limits -1 and +1. It is then possible to obtain (d Pg / dp) 
anal 5 d;ically from the polynomial coefficients. 

(4) Finally, the refractive index n(r) is obtained from equation (7,1) by 
numerical quadrature, Gaussian quadrature with 32 intervals is used 
after subdividing the range of integration into a number of subintervals. 

An example of a recovered lower atmospheric refractivity profile using the 
Herglotz-Wiechert algorithm is shown in Table III. The input profile is the same 
as that assumed for the first example in the iterative technique (Section 6.0). 
Recovery is .4% accurate at the surface and 40% accurate at 33 kilometers. 

Table III 

Results of Inversion by H-W Algorithm 
for the Three Parameter Profile 


Radius 

(km) 

Refractivity 
Computed from the 
H-W Algorithm 

Refractivity 
from the 
Input Profile 

Difference 

637^,02 

320.2 

321.4 

1.2 

6381.71 

212.2 

213.5 

1.3 

6384.13 

146.3 

147.7 

1.4 

6386.40 

103.1 

104.5 

1.4 

6388.59 

73.4 

74.9 

1.5 

6390.72 

52.7 

54.1 

1.4 

6392.82 

37.0 

39.3 

1.4 

6394.88 

27.3 

28.7 

1.4 

6396.93 

19.7 

21.0 

1.3 

6398.97 

14.1 

15.4 

1.3 

6401.00 

10.1 


1.2 

6403.01 

7.1 


1.2 

6405.03 

5.0 


1.1 

6407.04 

3.4 

4.5 

1.4 

6409.05 

2.3 

3.3 

1.0 

6411.05 

1.5 

2.4 

.9 
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8.0 RECOVERY OF PRESSURE AND TEMPERATURE 

Refractivity at radio frequencies has been previously defined (equation 3.1) as 


N = 77.6 (P/T) + (3.73)105 (e/T2) (8.1) 


If the wet component of N (second term on right side of (8.1)) can be properly 
modelled or determined (e.g. by an independent measurement), the remaining 
part, Nj^y, is proportional to density, a function of both pressure and tempera- 
ture 


Nj^y = 77.6 (P/T) (8.2) 


Assiuning hydrostatic equilibrimn and the perfect gas law, the pressure at a 
given height can be obtained by integrating Nj,y [ 21]. The temperature at the 
given height can then be calcvilated from (8.2). 

The hydrostatic equation is [ 13] 

dP = -p g(z) dz (8.3) 

where g(z) is ([22], pg. 217) 

g(z) = go 2 )^ (m/sec2) 

and gQ, the local acceleration of gravity is calculated from the latitude ([22], 
pg. 488) 

go = 9.780356 (1 + 0.0052885 sin^ " 5.9 x10-® sin ^ 2 (8.5) 
The effective earth's radius Tq in (8.4) is given by ([22], pg. 218) 

To = 2go/f(<^) (8.6) 
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where 


= 3. 085462 X 10“® + 2. 27 x lO'^ cos 2(^ - 2 x ^os 4 c^l 


In the above equations z is the geometric height and dz is the differential geo- 
metric height, both in geometric kilometers. 

The perfect gas law is ([ 13]) 

p = MP/RT (8,8) 

where 


M = the apparent molecular weight of dry air = 28.966 ([ 22], pg. 289) 

R = universal gas constant = 8314.36 Joules (®K)~^ x (Kg- mole)" ^ ([22], 
pg. 289) 

P = pressure in millibars 

= virtual temperature in degrees Kelvin 
P = density of dry air 
Substituting (8.8) in (8.3) and using (8.2) 



Taking the integral of both sides of (8.8), the pressure P at the height h is ob- 
tained 



M 


77.6 R 


J"g(z)(T/T,)Nj„ 
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or 


P 


(M/77.6R) 



g(z) (T/T^) (N 


NJdz 


( 8 . 10 ) 


The temperature T can be obtained using (8.2) and (8.10). 

From (8.10) it can be seen that, having recovered the total refractivity N, the 
wet component, must be subtracted out. In addition (T / T^) must be 

modelled. 

Using tables of the U.S. Standard Atmosphere Supplements [ 13] a regression 
model for both (T/ T^) and was obtained from 9 Standard atmospheres. 

For 15° North Annual and the 4 January atmospheres (30°N, 45°N, 60°N, 75°N), 
the following regression models were obtained 

In In (T^/T) = (-3.0557 - 0.0648 0 l) + [“0.6192 + 0.0043 h 

( 8 . 11 ) 

In \ = ( 5. 678 - 0. 0488 + [- 0. 5897 + (0. 1104) <^^3 h 


where is the latitude in degrees and h the altitude in kilometers. 

For the Jvily atmospheres (30°N, 45°N, 60°N, 75°N), the following regression 
models were obtained 


In In (T„/T) 
InN^ 


= (-3. 7414 - 0.0234 c^) - 0.4389 h 


( 8 . 12 ) 

= 5.4510 - 0.0186 0L" 0,4529 h - (0. 1756) lO'^ ^ 


Equation (8.10) can be written as 

P, = z) Nd(z) dz 

30 


(8.13) 



where 


K(0l. z) 




g(z) 


N 


d 


N-i 


and (T/ T^) and are given by either (8.11) or (8.12). 

Example 

From radiosonde data taken at Dulles airport in January 1967, the dry, wet, and 
total refractivity, N = were calculated from measured values of tem- 

perature, pressure and relative humidity using (3.1) and the procedure docu- 
mented in [ 15]. Both (T/ T„) and were estimated from (8.11) and the inte- 
grand in (8.10) munerically integrated from the height h to 200 kilometers to 
obtain P^. was obtained from (8.2). The results can be seen in Table IV. 
Recovery of pleasure is accurate to within 0.3%, The temperature recovery 
is worst at the surface, but improves with height. 
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Table IV 

Pressure and Temperature Recovery From Estimated Dry Refractivity 


Dry Calculated Pressure 

Height Refractivity Pressure 

(km) N^ = N-N„ P(mb) = | K(<^L-h)Nddh 100 

•'h 


270.67 1005 


254.53 

222.41 

199.19 
180.14 

163.19 
91.20 
19.94 


886.47 

787.69 

695.32 

609.28 

538.04 

265.45 

54.47 


1002.16 

885.95 

788.68 

696.84 

610.36 

538.51 

265.45 

54.47 


-0.29 

-0.06 

0.13 

0.22 

0.18 

0.09 

0.07 

0.27 


T (°K) 

77.6 

T = — * 

273.15 

287.32 

271.81 

270.10 

271.70 

275.18 

271.25 

271.48 

264.47 

262.93 

257.79 

256.08 

225.41 

225.86 

211.54 

211.55 



















9.0 SUMMARY 


Two techniques for recovering lower atmospheric refractivity profiles from 
simulated satellite-to-satellite tracking data have been described and examples 
given using the geometric configuration of the ATS-6/NIMBUS-6 Tracking Ex- 
periment. Both satellites are in circular orbit (NIMBUS in polar orbit and ATS 
in geostationary orbit in equatorial plane) and the underlying atmospheric refrac- 
tivity has the symmetric form N = exp P (s) where P (s) is a polynomial in the 
normalized height s. For purposes of recovery, a second degree polynomial was 
used for P(s). 

For the particular simulation used in the analysis, the Herglotz-Wiechert tech- 
nique recovered values which were 0.4% and 40% different from the input values 
(calculated from the parameters of the assumed model atmosphere) at the sur- 
face and at a height of 33 kilometers respectively. Using the same input data, the 
model fitting technique recovered refractivity values 0.05% and 1% different from 
the input values at the surface and at a height of 50 kilometers, respectively. 

It should be noted that modelling errors have not been included in the above. 

The examples cited show the relative degradation in the accuracy due to the 
calculation method rather than the true accuracy. 

The possible application of recovering atmospheric refractivity profiles lies in 
the area of meteorology. If the effects of the ionosphere and water vapor can be 
properly modelled or removed from the data, then calculation of pressure and 
temperature distributions can be made from the dry refractivity through numeri- 
cal integration procedures. This has been demonstrated using radiosonde data 
from Dulles airport, Virginia, calculating the total refractivity from measured 
values of temperature, pressure, and relative humidity, and using a regression 
model (obtained from 9 Standard atmospheres) to subtract out the wet refrac- 
tivity. The term (T/ T^) in the integrand of (8.10) was estimated from a regres- 
sional model (also obtained from 9 Standard atmospheres) and finally equation 
(8.10) numerically integrated from a height h to 200 kilometers to obtain pres- 
sure. The temperature was then obtained from (8.2). Pressure recovery was 
less than 0.3% to 20 kilometers. Temperature recovery at the surface was 5% 
but improved with height. 

The techniques documented here will be applied to the calculation of lower atmo- 
spheric refractivity profiles from actual satellite-to-satellite tracking data taken 
during the radio occultation portion of the ATS-6/ NIMBUS-6 Tracking and Orbit 
Determination Experiment. The basic approach will be the iterative model fit- 
ting technique, using the Herglotz-Wiechert technique as a "starter" to obtain 
an initial estimate of the parameter vector (apriori values). It is planned to ob- 
tain radiosonde data as close as possible to the point viiefe the ray path of the 
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signal between the two satellites comes closest to the earth at the time of occul- 
tation. From the radiosonde data refractivity profiles will be calculated to 
verify and compare against those obtained by radio occultation methods. In addi- 
tion both pressure and temperature distributions will be obtained using the pro- 
cedure described above and compared against the radiosonde data. 

Along with the analysis of the occultation data from ATS/NIMBUS there will be 
a continuing effort to model the water vapor and investigate horizontal gradients 
using regression techniques. 


10. CONCLUSIONS 

(1) Two techniques for recovering lower atmospheric refractivity profiles from 
simulated satellite-to- satellite tracking data have been documented — the 
Herglotz-Wiechert algorithm employing the integral Abel transform and an 
iterative model fitting method incorporating a numerical ray tracing proce- 
dure. Of the two the model fitting technique recovered the input profile 
more accurately. 

(2) For accurate inversion of lower atmospheric refractivity profiles, it is 
necessary to take the ionospheric refraction of the ray into account. A 
simple model of the ionospheric electron density like the modified Chapman 
profile has been shown to be adequate for accurate recovery of the lower 
atmospheric profile. 

(3) By properly modelling the ionosphere and water vapor effects (through re- 
gression techniques), it is possible to obtain pressure and temperature 
distributions from the dry refractivity using the hydrostatic equation and 
the perfect gas law. 
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APPENDIX A 


WEIGHTED NONLINEAR LEAST SQUARES WITH APRIORI 
PARAMETER VECTOR AND COVARIANCE MATRIX^ 


The particular numerical method used to solve the nonlinear least squares 
problem is a modified Newton-Raphson iteration technique. 

Let r be an (nx 1) vector of measurements andf(x) a nonlinear function or 
model v/hich it is assumed adequately describes the measurements, where x is 
an (m x 1) parameter vector. 


( Z J, Z2, . . 

•> ^n) 

(Xi, Xj, . .. 



(A.1) 


It should be noted that in the occxdtation problem where the refractivity has the 
form N - exp P(s) and P(s) = ag + a^ s + ajS^ + . . . + a^_^s"’"^ is a polynomial 
in the normalized height s = (.01) h - 1 (h is height in kilometers), x is the 
(m X 1) parameter vector of coefficients a^. The vector z is the vector of ob- 
served ranges between the two satellites or the vector of observed range rates 
along the propagation path of the signal, or both of these quantities. The function 
f (x) is a function of the parameter vector x which can be the calculated range , 
range rate, or both. In any case, the parameters a. enter in a nonlinear fashion 
in the calculation of these quantities (ray tracing of angles and distances within 
the troposphere). We can write 

= (fi- fa fn) (A.2) 


where f . represents the value of range (or range rate) calculated for the i^ data 
point for a particular set of values for x (i.e. , the a /s). 

^ This appendix is self contained. The notation is therefore general. 
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Let X be an apriori estimate of the parameter vector x with its associated 
covariance matrix Dj^ (m x m matrix). Let D, be the covariance matrix of the 
observation vector z . The function that we wish to minimize is the quadratic 
form 


Q - [z - f (x [z - f(x)] + (x -x)'*‘DgT*(x -x) (A.3) 

where 

X - (Xj, Xj, •• • . (A. 4) 

In a statistical sense we can consider the apriori values as random variables 
(or additional observations) having the mean vector x. The covariance matrix 
then represents the confidence or certainty that we have in our "best" esti- 
mate X of the parameter vector x . As an example of how one might arrive at 
values for x and suppose that from tabular values of pressure, temperature, 
and relative humidity (either radiosonde data from a particular weather station 
or data from a standard atmosphere) a refractivity profile is calculated as a 
fimction of height. Using linear least squares we can then fit the logarithm of 
the calculated values of refractivity with a polynomial of degree (m - 1), The 
solution vector of coefficients a^ and its covariance matrix could then be used 
as an apriori parameter vector "x with covariance matrix Dj^ for a subsequent 
occiiltation problem. In another sense the matrix serves as a control in the 
iteration process. From (A.3) it can be seen that if Djj. is very large (large im- 
certainty in x ), then will be very small and the second term on the right 
side of (A.3) becomes negligible, which means that the observations "take over" 
and determine the solution (standard least squares problem); if the opposite is 
true, i.e., is very small, the solution will be near the apriori vector x. 

Thus, serves as a weighting matrix for the apriori values. In a similar 
manner we can talk about the matrix which is usually taken as diagonal (ob- 
servations are uncorrelated). For least squares with equal weighting is a 
scalar times the identity matrix I. 

In order for (A.3) to be minimized, it is necessary that the gradient of Q with 
respect to X , Q, evaluated at the solution X = x*, vanish 

F(x) = V^Q = [z - f (x)] -2D/1 (x - x) (A.5) 
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where 



B(x) 


(A.6) 


is an (n X m) matrix of partial derivatives of the calculated quantities with re- 
spect to the parameters of rank m < n. 

Neglectii^ the dependence of B on x a modified Newton-Raphson iteration tech- 
nique can be used to find a solution x* such that F (x *) = 0, Taking the gradient 
of F(x) with respect to x 


F'(x) - + 2D5f ^ 


(A.7) 


Since both Dj and are covariance matrices, they are by definition positive 
definite. Since B is assumed to be of rank m < n, B^ D^'^B is positive definite. 
Therefore (A.7) is positive definite, and if the process converges to x*, Q will 
be minimized. 

In order to solve for x*, denote by x^"*^ the estimate of the solution at the m^ 
iteration; then, the improved estimate of the solution x = x’*‘, using the 
Newton-Raphson formula will be given by 

X* = x<”*^> = x("*> - F(x('"))/F'(x('")) 


or 

X* = x^") + (B^D^^B + jB^D-^ [z - f (x W)] + x« (A.8) 

where 

B - B(x<"’>) 
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Assuming l^e observations z to be independent (in the statistical sense) of the 
apriori values x , the covariance matrix of the solution x’X is given by 

+ DjtI)-* (A.9) 

The Uncertainty in Refractivity 
Let the refractivity model be given by 

N = ejq) s^x* = e:qp P(s) (A.10) 


where x* is the solution vector 


and 


(x*)’r 




- (1 S ... S™"*) (■^•11) 

is an (m X 1) vector of normalized heights s = (O.Ol)h - 1 (h is height in 
kilometers). 

An error AN in the neighborhood of the solution x * is 

an = exp[s^(x* + Ax)] - (s‘*'x*) 

= [exp (s’^'x*)] s^Ax (A.12) 


The variance of AN (or of N) is, using (A.9) 

cr^j2 = exp (2s’‘x*) [b^Dj-^ + DgT^] ^ s (A. 13) 

The uncertainty or standard deviation of N is the square root of (A.13). 
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APPENDIX B 


RAY TRACING THROUGH THE LOWER ATMOSPHERE AND IONOSPHERE 


From Fermat's principle of extremal path and the Euler- Lagrange variational 
equation, we obtain for a ray, 


n(r) • r sin i - p 


(B.l) 


where n is the phase refractive index, r the geocentric radius, i the angle of 
incidence and p the distance of the ray asymptote from the earth's center. The 
total bending angle r (Figure 10) is given by 


T = -2 



dn(r) dr 



(B. 2 ) 


where 77 = n (r) • r 

r^ = radius of closest approach of the ray 
r^ = radius of the outer limit of the atmosphere 

The angle subtended by the total arc of the ray in the atmosphere at the center 
of the earth 4 > is given by 


4 > = r + 77- - 2i^ 


(B. 3 ) 


where i^ is the angle of incidence at the outer limit of the atmosphere given by 

i. = 


The computed value of the angle subtended by the total ray path at the center of 
the earth is given by 



4>^ cj> 


(B.5) 


where 4>^ and are the angles subtended at the center of the earth by the linear 
segments of the ray. 

The observed angle subtended by the total 
given by 


(9q - cos"^ 


where and R^ are the radius vectors of NIMBUS and ATS respectively. The 
ray linking NIMBUS and ATS for a given atmosphere is found by iteratively 
altering rp until d ^ agrees with 8^ to within the required precision. 

After finding the specific ray linking ATS and NIMBUS, the one-way raiige along 
that ray is computed from 


ray path at the center of the earth is 

Rn • «A 

^ Ra 


¥ 

(B.6) 


Rp = NB + CA + P (B.7) 

where NB and CA are the linear segments of the ray which can be calculated 
from a knowledge of Rj^, R^, r^, and i„. P is the path length of the ray in the 
atmosphere. 

The lower atmosphere is non-dispersive for radio waves and hence P is given by 



t 
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The ionosphere is dispersive and here P is given by 


P 


I 



r>g (r) dr 
cos i 


where n g is the group refractive index. 


(B.9) 




and f is the frequency of transmission, n is the phase index of refraction. At 
S-band frequencies (2 MHz), P can be written as 



(B.11) 


where r j is the radius of the upper limit of the ionosphere (altitude of 1000 
kilometers for this analysis). 

The integrals (B.2), (B.9), and (B.ll) are evaluated by Gaussian quadrature. 

The range rate between the two satellites is obtained by the algebraic addition 
of the component of the velocity of each satellite along the ray path at the satel- 
lite. The ray path direction at each satellite is found by knowing the angle be- 
tween the ray path and the radial line from the center of the earth to the satellite. 
For NIMBUS (Figure 10) this angle '/'f, is 


= sin-^ l~-\ (B.^l^^^^^ 

\ j 
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Similarly, the angle between the ray path and the radius at ATS t/>y^ is given by 



(B.13) 



APPENDIX C 

HERGLOTZ-WIECHERT ALGORITHM FOR INVERSION OF PHASE DATA 


From Fermat's principle of extremal path and the Euler- Lagrange variational 
equation we obtain for a ray 


t 


n(r) • r sin i - p 


(C.l) 


where p is the distance of the ray asymptote from the earth's center, n is the 
phase index of refraction at a radial distance r and i is the angle of incidence. 

The phase path length P along the ray is given by 


P(yo) = 


f ” ' ^ (— ) d7? 

K r ^772 - p2 VdTj/ 


where 77 = n(r) • r 

and r^ is tiie radius to the outer limit of the atmosphere. 
Transforming variables 



we obtain 


(C.2) 


W ( 2 _ y) [d (In r)/dy] 

P(w) = -2 I dy 

*0 yw - 7 


(G.3) 


G-1 



The Abel transform of P(w) is 


.•(y) = (r.^ - y> ^ 


(C.4) 


or 


g'(y) " 


dy 


kj: 


f(w) 




dw 

w 


(C.5) 


Transforming back from y and w to rj and p respectively 


d(lnr) 

1 ^ 

1 

r’ pP(p) dp 

dr} 

1 7^2 1 dr] 

rr 



(C.6) 


By integrating (C.6), changing the order of integration and by successive integra- 
tion by parts we obtain 


r('^) = r^exp 


' I r’ 

J oosh-^ ip/v) (1//0) 


dP dp 
dp 


(C.7) 


n(r) - T?/r(77) 


(C.8) 


At the turning point of the ray 


sin i - 1 
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p 


(C.9) 


\ = 


p p 


where the subscript p refers to the turning point. 
By combining (C.7), (C. 8 ), and (C.9) we get 


"p = SI" Im 


I 


(1/tt) I cosh“* (p/t?) — 


1 dPdp 


P Ap 


(C.IO) 


By considering the total phase change as the sum of the phase change due to 
the relative motion of the satellites in vacuum and the phase change P 3 due to 
the atmosphere we can write 


P(P) = P + P 3 


(C.ll) 


Hence, in the absence of an atmosphere 


P(P) = Pe 


and 


sin exp 


r 

•'r 


1 dP dp 

(I/77) I co^-^ (p/t))— — -- 


(C.12) 


Hence, we can rewrite (C.IO) as 


ejqp 


(1/77) 




1 dp3 dp 

7 -77- 


(C.13) 
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